

###########################################################################
###########################################################################
## ZIMBABWE ANALYSIS - ANGER FIELD EXPERIMENT
## 2021 Apr 01
## Lauren Young
###########################################################################
###########################################################################

rm(list=ls())

library(data.table)
library(stargazer)
library(gtools)
library(lubridate)
library(foreign)
library(xtable)
require(ggplot2)  
library(car)
library(interflex)
library(maptools)
library(ri2)
library(sandwich)
library(lmtest)
library(ggspatial)
library(cowplot)


## set working directory

setwd('~/Box/Zimbabwe/party_experiment/04_Analysis/replication_files')

## set seed

set.seed(103648)


## read in functions

st <- function(x) {
  (x - mean(x,na.rm=T))/sd(x,na.rm=T)
}

st_out <- function(x, data = dat, treat = 'treat_all') {
  (x - mean(x[data[,treat]==0],na.rm=T))/sd(x[data[,treat]==0],na.rm=T)
}


## read in data

dat <- read.dta("01_data/02_clean/dat_message.dta")
dat2 <- read.csv("01_data/02_clean/dat_group.csv")



##########
## SUMMARY STATS
##########

#####
## text - compliance stats (p20)
#####

table(dat2$t1, dat2$x)
table(dat2$prov, dat2$x)


#####
## Figure A1 - response over time
#####

tab <- data.table(dat)
tab <- tab[,list('msg_count' = sum(message), 
                 'viva_count' = sum(viva), 'vic_count' = sum(vic),
                 'resp_count' = length(unique(id))),
                      by = list(t1, admin, treat_dur_hr)]
tab <- tab[tab$admin==0 & is.na(tab$t1)==F,]

par(mfrow=c(1,1))

pdf('03_graphs/message_time.pdf')
par(mar=c(4,6,4,2))
plot(y=tab$msg_count[tab$t1==1], x=tab$treat_dur[tab$t1==1],
     type = 'l', col = 'red',
     xlab = 'Hours from Treatment', ylab = 'Messages',
     cex.axis = 2.4, cex.lab = 2.4, lwd = 4, lty = 1)
points(y=tab$msg_count[tab$t1==0], x=tab$treat_dur[tab$t1==0],
       type = 'l', col = 'darkgrey', lty = 6, lwd = 4)
abline(v=0, lty = 2)
legend('topleft', c('anger', 'enthusiasm'), col = c('red', 'darkgrey'), lty = c(1,6), lwd = 4, cex = 2)
dev.off()

pdf('03_graphs/resp_time.pdf')
par(mar=c(4,6,4,2))
plot(y=tab$resp_count[tab$t1==1], x=tab$treat_dur[tab$t1==1],
     type = 'l', col = 'red',
     xlab = 'Hours from Treatment', ylab = 'Respondents',
     cex.axis = 2.4, cex.lab = 2.4, lwd = 4, lty = 1)
points(y=tab$resp_count[tab$t1==0], x=tab$treat_dur[tab$t1==0],
       type = 'l', col = 'darkgrey', lty = 6, lwd = 4)
abline(v=0, lty = 2)
legend('topleft', c('anger', 'enthusiasm'), col = c('red', 'darkgrey'), lty = c(1,6), lwd = 4, cex = 2)
dev.off()

pdf('03_graphs/viva_time.pdf')
par(mar=c(4,6,4,2))
plot(y=tab$viva_count[tab$t1==1], x=tab$treat_dur[tab$t1==1],
     type = 'l', col = 'red',
     xlab = 'Hours from Treatment', ylab = 'Slogans',
     cex.axis = 2.4, cex.lab = 2.4, lwd = 4, lty = 1)
points(y=tab$viva_count[tab$t1==0], x=tab$treat_dur[tab$t1==0],
       type = 'l', col = 'darkgrey', lty = 6, lwd = 4)
abline(v=0, lty = 2)
legend('topleft', c('anger', 'enthusiasm'), col = c('red', 'darkgrey'), lty = c(1,6), lwd = 4, cex = 2)
dev.off()

pdf('03_graphs/vic_time.pdf')
par(mar=c(4,6,4,2))
plot(y=tab$vic_count[tab$t1==1], x=tab$treat_dur[tab$t1==1],
     type = 'l', col = 'red',
     xlab = 'Hours from Treatment', ylab = 'Symbols',
     cex.axis = 2.4, cex.lab = 2.4, lwd = 4, lty = 1)
points(y=tab$vic_count[tab$t1==0], x=tab$treat_dur[tab$t1==0],
       type = 'l', col = 'darkgrey', lty = 6, lwd = 4)
abline(v=0, lty = 2)
legend('topleft', c('anger', 'enthusiasm'), col = c('red', 'darkgrey'), lty = c(1,6), lwd = 4, cex = 2)
dev.off()


#####
## text - count total response by treatment assignment (p27)
#####

dat <- data.table(dat)
tab <- dat[,list('msg_count' = sum(message), 
                 'viva_count' = sum(viva), 'vic_count' = sum(vic),
                 'resp_count' = length(unique(id))),
           by = list(t1, admin)]
tab <- tab[tab$admin==0 & is.na(tab$t1)==F,]
tab


#####
## TABLE A1 - balance
#####

tab <- data.frame('Variable' = c('Pre-Treatment Messages', 'Pre-Treatment Respondents',
                            'Pre-Treatment Slogans', 'Pre-Treatment Symbols',
                            'Poverty', 'Violence - Sokwanele',
                            'Violent Events - ACLED', 'Fatalities - ACLED'),
                  'Control Mean' = c(mean(dat2$msgs.pre[dat2$t1==0]),
                                       mean(dat2$resp.pre[dat2$t1==0]),
                                       mean(dat2$viva.pre[dat2$t1==0]),
                                       mean(dat2$vic.pre[dat2$t1==0]),
                                       mean(dat2$wgt_z[dat2$t1==0]),
                                       mean(dat2$any_bin[dat2$t1==0]),
                                     mean(dat2$acled_events_0015[dat2$t1==0]),
                                     mean(dat2$acled_fatalities_0015[dat2$t1==0])),
                  'Treatment Mean' = c(mean(dat2$msgs.pre[dat2$t1==1]),
                                     mean(dat2$resp.pre[dat2$t1==1]),
                                     mean(dat2$viva.pre[dat2$t1==1]),
                                     mean(dat2$vic.pre[dat2$t1==1]),
                                     mean(dat2$wgt_z[dat2$t1==1]),
                                     mean(dat2$any_bin[dat2$t1==1]),
                                     mean(dat2$acled_events_0015[dat2$t1==1]),
                                     mean(dat2$acled_fatalities_0015[dat2$t1==1])),
                  'p-value' = c(t.test(dat2$msgs.pre[dat2$t1==0], dat2$msgs.pre[dat2$t1==1])$p.value,
                                t.test(dat2$resp.pre[dat2$t1==0], dat2$resp.pre[dat2$t1==1])$p.value,
                                t.test(dat2$viva.pre[dat2$t1==0], dat2$viva.pre[dat2$t1==1])$p.value,
                                t.test(dat2$vic.pre[dat2$t1==0], dat2$vic.pre[dat2$t1==1])$p.value,
                                t.test(dat2$wgt_z[dat2$t1==0], dat2$wgt_z[dat2$t1==1])$p.value,
                                t.test(dat2$any_bin[dat2$t1==0], dat2$any_bin[dat2$t1==1])$p.value,
                                t.test(dat2$acled_events_0015[dat2$t1==0], dat2$acled_events_0015[dat2$t1==1])$p.value,
                                t.test(dat2$acled_fatalities_0015[dat2$t1==0], dat2$acled_fatalities_0015[dat2$t1==1])$p.value))
tab
rownames(tab) <- tab$Variable; tab$Variable=NULL
xtable(tab)




#####
## test of intercoder reliability
#####

## anger
t <- table(dat$anger_c1, dat$anger_c2)
t
(t['No','No'] + t["Yes","Yes"])/(t['No','No'] + t["Yes","Yes"] + t['No','Yes'] + t["Yes","No"])

## enthusiasm
t <- table(dat$enthusiasm_c1, dat$enthusiasm_c2)
t
(t['No','No'] + t["Yes","Yes"])/(t['No','No'] + t["Yes","Yes"] + t['No','Yes'] + t["Yes","No"])
prop.table(table(dat$enthusiasm_c1))
prop.table(table(dat$enthusiasm_c2))

## sadness
t <- table(dat$sadness_c1, dat$sadness_c2)
t
(t['No','No'] + t["Yes","Yes"])/(t['No','No'] + t["Yes","Yes"] + t['No','Yes'] + t["Yes","No"])

## fear
t <- table(dat$fear_c1, dat$fear_c2)
t
(t['No','No'] + t["Yes","Yes"])/(t['No','No'] + t["Yes","Yes"] + t['No','Yes'] + t["Yes","No"])




#####
## Table 1 - manipulation check 
#####


## ITTs

mod0a <- lm(anger_any ~ t1 + block, data = dat2, weights = ipw)
mod0a$se <- coeftest(mod0a, vcov = vcovHC, type = "HC1")[,2]
mod0b <- lm(anger_any ~ t1 + block + members_log + any_bin + wgt_z, data = dat2, weights = ipw)
mod0b$se <- coeftest(mod0b, vcov = vcovHC, type = "HC1")[,2]
mod1a <- lm(enth_any ~ t1 + block, data = dat2, weights = ipw)
mod1a$se <- coeftest(mod1a, vcov = vcovHC, type = "HC1")[,2]
mod1b <- lm(enth_any ~ t1 + block + members_log + any_bin + wgt_z, data = dat2, weights = ipw)
mod1b$se <- coeftest(mod1b, vcov = vcovHC, type = "HC1")[,2]
mod2a <- lm(sadness_any ~ t1 + block, data = dat2, weights = ipw)
mod2a$se <- coeftest(mod2a, vcov = vcovHC, type = "HC1")[,2]
mod2b <- lm(sadness_any ~ t1 + block + members_log + any_bin + wgt_z, data = dat2, weights = ipw)
mod2b$se <- coeftest(mod2b, vcov = vcovHC, type = "HC1")[,2]
mod3a <- lm(fear_any ~ t1 + block, data = dat2, weights = ipw)
mod3a$se <- coeftest(mod3a, vcov = vcovHC, type = "HC1")[,2]
mod3b <- lm(fear_any ~ t1 + block + members_log + any_bin + wgt_z, data = dat2, weights = ipw)
mod3b$se <- coeftest(mod3b, vcov = vcovHC, type = "HC1")[,2]

stargazer(mod0a, mod0b, mod1a, mod1b, mod2a, mod2b, mod3a, mod3b, 
          se = list(mod0a$se, mod0b$se, mod1a$se, mod1b$se, mod2a$se, mod2b$se, mod3a$se, mod3b$se),
          digits = 2, no.space = T, keep.stat = c('n', 'rsq'), 
          covariate.labels = c('Anger Appeal', 'Members (Log)', 
                               'Any Violence', 'Poverty', 
                               'Constant'),
          omit = c('block'))



#####
## Table 2 - outcomes 
#####

## ITTs

dat2$pre <- dat2$index.pre.m
mod0a <- lm(index.m ~ t1 + block, data = dat2, weights = ipw)
mod0a$se <- coeftest(mod0a, vcov = vcovHC, type = "HC1")[,2]
mod0b <- lm(index.m ~ t1 + block + pre + members_log + any_bin + wgt_z, data = dat2, weights = ipw)
mod0b$se <- coeftest(mod0b, vcov = vcovHC, type = "HC1")[,2]
dat2$pre <- dat2$msgs.pre
mod1a <- lm(msgs ~ t1 + block, data = dat2, weights = ipw)
mod1a$se <- coeftest(mod1a, vcov = vcovHC, type = "HC1")[,2]
mod1b <- lm(msgs ~ t1 + block + pre + members_log + any_bin + wgt_z, data = dat2, weights = ipw)
mod1b$se <- coeftest(mod1b, vcov = vcovHC, type = "HC1")[,2]
dat2$pre <- dat2$resp.pre
mod2a <- lm(resp ~ t1 + block, data = dat2, weights = ipw)
mod2a$se <- coeftest(mod2a, vcov = vcovHC, type = "HC1")[,2]
mod2b <- lm(resp ~ t1 + block + pre + members_log + any_bin + wgt_z, data = dat2, weights = ipw)
mod2b$se <- coeftest(mod2b, vcov = vcovHC, type = "HC1")[,2]
dat2$pre <- dat2$viva.pre
mod3a <- lm(viva ~ t1 + block, data = dat2, weights = ipw)
mod3a$se <- coeftest(mod3a, vcov = vcovHC, type = "HC1")[,2]
mod3b <- lm(viva ~ t1 + block + pre + members_log + any_bin + wgt_z, data = dat2, weights = ipw)
mod3b$se <- coeftest(mod3b, vcov = vcovHC, type = "HC1")[,2]
dat2$pre <- dat2$vic.pre
mod4a <- lm(vic ~ t1 + block, data = dat2, weights = ipw)
mod4a$se <- coeftest(mod4a, vcov = vcovHC, type = "HC1")[,2]
mod4b <- lm(vic ~ t1 + block + pre + members_log + any_bin + wgt_z, data = dat2, weights = ipw)
mod4b$se <- coeftest(mod4b, vcov = vcovHC, type = "HC1")[,2]

stargazer(mod0a, mod0b, mod1a, mod1b, mod2a, mod2b, mod3a, mod3b, mod4a, mod4b,
          se = list(mod0a$se, mod0b$se, mod1a$se, mod1b$se, mod2a$se, mod2b$se, mod3a$se, mod3b$se, mod4a$se, mod4b$se),
          digits = 2, no.space = T, keep.stat = c('n', 'rsq'), 
          covariate.labels = c('Anger Appeal', 'Pre-Treatment Outcome',
                               'Members (Log)', 
                               'Any Violence', 'Poverty', 
                               'Constant'),
          omit = c('block'))



#####
## Text Section 5.2 - percent change in treated
#####

ates <- c(mod1b$coef['t1'],mod2b$coef['t1'],mod3b$coef['t1'],mod4b$coef['t1'])
means1 <- c(mean(dat2$msgs[dat2$t1==1],na.rm=T), mean(dat2$resp[dat2$t1==1],na.rm=T),
            mean(dat2$viva[dat2$t1==1],na.rm=T), mean(dat2$vic[dat2$t1==1],na.rm=T))
means0 <- c(mean(dat2$msgs[dat2$t1==0],na.rm=T), mean(dat2$resp[dat2$t1==0],na.rm=T),
            mean(dat2$viva[dat2$t1==0],na.rm=T), mean(dat2$vic[dat2$t1==0],na.rm=T))

tab <- data.frame('means0' = means0, 'means1' = means1, 'diff' = ates/means0)
rownames(tab) <- c('Messages', 'Respondents', 'Slogans', 'Symbols')
xtable(tab)



#####
## Figure 4 - ATEs over time
#####

## make an index for each 3-hour window

for (i in seq(3,21,3)){
  i2 <- (dat2[,paste0('msgs',i)] - mean(dat2[,paste0('msgs',i)], na.rm=T))/sd(dat2[,paste0('msgs',i)][dat2$t1==0], na.rm=T)
  i3 <- (dat2[,paste0('resp',i)] - mean(dat2[,paste0('resp',i)], na.rm=T))/sd(dat2[,paste0('resp',i)][dat2$t1==0], na.rm=T)
  i4 <- (dat2[,paste0('viva',i)] - mean(dat2[,paste0('viva',i)], na.rm=T))/sd(dat2[,paste0('viva',i)][dat2$t1==0], na.rm=T)
  i5 <- (dat2[,paste0('vic',i)] - mean(dat2[,paste0('vic',i)], na.rm=T))/sd(dat2[,paste0('vic',i)][dat2$t1==0], na.rm=T)
  dat2$temp <- (i2 + i3 + i4 + i5)/5
  setnames(dat2, 'temp', paste0('index',i))
  i2 <- (dat2[,paste0('msgs',i,'.pre')] - mean(dat2[,paste0('msgs',i,'.pre')], na.rm=T))/sd(dat2[,paste0('msgs',i,'.pre')][dat2$t1==0], na.rm=T)
  i3 <- (dat2[,paste0('resp',i,'.pre')] - mean(dat2[,paste0('resp',i,'.pre')], na.rm=T))/sd(dat2[,paste0('resp',i,'.pre')][dat2$t1==0], na.rm=T)
  i4 <- (dat2[,paste0('viva',i,'.pre')] - mean(dat2[,paste0('viva',i,'.pre')], na.rm=T))/sd(dat2[,paste0('viva',i,'.pre')][dat2$t1==0], na.rm=T)
  i5 <- (dat2[,paste0('vic',i,'.pre')] - mean(dat2[,paste0('vic',i,'.pre')], na.rm=T))/sd(dat2[,paste0('vic',i,'.pre')][dat2$t1==0], na.rm=T)
  i2 <- ifelse(i2==Inf | i2==-Inf, NA, i2)
  i3 <- ifelse(i3==Inf | i3==-Inf, NA, i3)
  i4 <- ifelse(i4==Inf | i4==-Inf, NA, i4) 
  i5 <- ifelse(i5==Inf|-Inf, NA, i5)
  dat2$temp <- rowMeans(cbind(i2,i3,i4,i5),na.rm=T)
  setnames(dat2, 'temp', paste0('index',i,'.pre'))
}


## calculate ITTs

dat2$index24 <- dat2$index.m
outs = c(paste0('index', seq(3,24,3)))
treats = c('t1')

tab <- expand.grid('outs' = outs, 'treats' = treats)
tab$outs <- as.character(tab$outs)
tab$ate <- NA
tab$p <- NA
tab$se_lo <- NA
tab$se_hi <- NA


for (i in 1:dim(tab)[1]){
  dat2$outcome <- dat2[,tab$outs[i]]
  mod <- lm(as.formula(paste0(tab$outs[i], ' ~ t1 + block + members_log + any_bin + wgt_z')),
                       data = dat2, weights = ipw)
  tab$ate[i] <- summary(mod)$coefficients['t1',1]
  tab$p[i] <- summary(mod)$coefficients['t1',4]
  tab$se_lo[i] <- summary(mod)$coefficients['t1',1]-1.96*summary(mod)$coefficients['t1',2]
  tab$se_hi[i] <- summary(mod)$coefficients['t1',1]+1.96*summary(mod)$coefficients['t1',2]
}

tab$period <- seq(3, 24, 3)

pdf('03_graphs/ate_time_index.pdf')
ggplot(tab, aes(period, ate)) +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
  geom_errorbar(aes(x = period, y = ate, ymin = se_lo, ymax = se_hi),
                lwd = 1, position = position_dodge(width = 1/2),
                width = 0.35, colour = "red") +
  geom_point(aes(x = period, y = ate), 
             position = position_dodge(width=1/2), colour = "red") +
  theme_minimal() +
  theme(text = element_text(size = 20),
        axis.text.x = element_text(angle = 0, hjust = 0.5)) +
  scale_y_continuous(name = 'Standardized Coefficient') + 
  scale_x_continuous(name = "Hours Post-Treatment", seq(3, 24, 3))
dev.off()



#####
## Table A3 (appendix) - Heterogeneous effects
#####

dat2$members_log_st <- st(dat2$members_log); dat2$membersXt1 <- dat2$members_log_st * dat2$t1

## models

dat2$pre <- st(dat2$index.pre.m)
mod0a <- lm(index.m ~ t1 + t2 + t1:t2 + t1:poor + t1:any_bin + any_bin + poor + block,
            data = dat2, weights = ipw)
mod0a$se <- coeftest(mod0a, vcov = vcovHC, type = "HC1")[,2]
mod0b <- lm(index.m ~ t1 + t2 + t1:t2 + t1:poor + t1:any_bin + any_bin + poor + block
            + members_log_st + pre, data = dat2, weights = ipw)
mod0b$se <- coeftest(mod0b, vcov = vcovHC, type = "HC1")[,2]
dat2$pre <- dat2$msgs.pre
mod1a <- lm(msgs ~ t1 + t2 + t1:t2 + t1:poor + t1:any_bin + any_bin + poor + block,
            data = dat2, weights = ipw)
mod1a$se <- coeftest(mod1a, vcov = vcovHC, type = "HC1")[,2]
mod1b <- lm(msgs ~ t1 + t2 + t1:t2 + t1:poor + t1:any_bin + any_bin + poor + block
            + members_log_st + pre, data = dat2, weights = ipw)
mod1b$se <- coeftest(mod1b, vcov = vcovHC, type = "HC1")[,2]
dat2$pre <- dat2$resp.pre
mod2a <- lm(resp ~ t1 + t2 + t1:t2 + t1:poor + t1:any_bin + any_bin + poor + block,
            data = dat2, weights = ipw)
mod2a$se <- coeftest(mod2a, vcov = vcovHC, type = "HC1")[,2]
mod2b <- lm(resp ~ t1 + t2 + t1:t2 + t1:poor + t1:any_bin + any_bin + poor + block
            + members_log_st + pre, data = dat2, weights = ipw)
mod2b$se <- coeftest(mod2b, vcov = vcovHC, type = "HC1")[,2]
dat2$pre <- dat2$viva.pre
mod3a <- lm(viva ~ t1 + t2 + t1:t2 + t1:poor + t1:any_bin + any_bin + poor + block,
            data = dat2, weights = ipw)
mod3a$se <- coeftest(mod3a, vcov = vcovHC, type = "HC1")[,2]
mod3b <- lm(viva ~ t1 + t2 + t1:t2 + t1:poor + t1:any_bin + any_bin + poor + block
            + members_log_st + pre, data = dat2, weights = ipw)
mod3b$se <- coeftest(mod3b, vcov = vcovHC, type = "HC1")[,2]
dat2$pre <- dat2$vic.pre
mod4a <- lm(vic ~ t1 + t2 + t1:t2 + t1:poor + t1:any_bin + any_bin + poor + block,
            data = dat2, weights = ipw)
mod4a$se <- coeftest(mod4a, vcov = vcovHC, type = "HC1")[,2]
mod4b <- lm(vic ~ t1 + t2 + t1:t2 + t1:poor + t1:any_bin + any_bin + poor + block
            + members_log_st + pre, data = dat2, weights = ipw)
mod4b$se <- coeftest(mod4b, vcov = vcovHC, type = "HC1")[,2]

stargazer(mod0a, mod0b, mod1a, mod1b, mod2a, mod2b, mod3a, mod3b, mod4a, mod4b,
          se = list(mod0a$se, mod0b$se, mod1a$se, mod1b$se, mod2a$se, mod2b$se, mod3a$se, mod3b$se, mod4a$se, mod4b$se),
          omit = c("block"),
          covariate.labels = c('Anger Appeal', 'Power Treatment', 'Any Violence',
                               'Poverty', 'Members (Log)', 'Pre-Treat Outcome',
                               'Anger Appeal $\\times$ Power Treatment',
                               'Anger Appeal $\\times$ Poverty', 
                               'Anger Appeal $\\times$ Any Violence', 
                               'Constant'),
          no.space = T, digits = 2,
          keep.stat = c('n', 'rsq'))



#####
## Figure 5 - Marginal effects plots
#####

dat2$viol <- dat2$any_bin; dat2$violXt1 <- dat2$viol * dat2$t1
dat2$poorXt1 <- dat2$poor * dat2$t1
dat2$members_log_st <- st(dat2$members_log); dat2$membersXt1 <- dat2$members_log_st * dat2$t1
dat2$t1Xt2 <- dat2$t1 * dat2$t2

dat2$pre <- dat2$index.pre.m

pdf('03_graphs/marg_viol_index.pdf')
print(inter.binning(Y = 'index.m', D = "t1", X = 'viol',
                    data = dat2,  cutoffs = c(0.5, 1.5),
                    Z = c('poor', 't2', 't1Xt2', 'poorXt1', 'pre', 'block'),
                    weights = 'ipw',
                    na.rm = T, vartype = 'robust', figure = T,
                    Ylabel = "Mean Effects Index", Dlabel = "Anger Treatment", Xlabel="Violence",
                    ylab = "Marginal Effect of Anger Treatment",
                    show.grid = T, theme.bw = T, cex.lab = 1.5, cex.axis = 1.5))
dev.off()



pdf('03_graphs/marg_poor_index.pdf')
print(inter.binning(Y = 'index.m', D = "t1", X = 'poor',
                    data = dat2, cutoffs = c(0,1),
                    Z = c('viol', 't2', 't1Xt2', 'violXt1', 'pre', 'block'),
                    weights = 'ipw', 
                    na.rm = T, vartype = 'robust', figure = T,
                    Ylabel = "Mean Effects Index", Dlabel = "Anger Treatment", Xlabel="Poverty",
                    ylab = "Marginal Effect of Anger Treatment",
                    show.grid = T, theme.bw = T, cex.lab = 1.5, cex.axis = 1.5))
dev.off()

pdf('03_graphs/marg_t2_index.pdf')
print(inter.binning(Y = 'index.m', D = "t1", X = 't2',
                    data = dat2, cutoffs = c(0.5, 1.5),
                    Z = c('viol', 'violXt1', 'poor', 'poorXt1', 'pre', 'block'), 
                    weights = 'ipw', 
                    na.rm = T, vartype = 'robust', figure = T,
                    Ylabel = "Mean Effects Index", Dlabel = "Anger Treatment", Xlabel="Power Treatment",
                    ylab = "Marginal Effect of Anger Treatment",
                    show.grid = T, theme.bw = T, cex.lab = 1.5, cex.axis = 1.5))
dev.off()



#####
## Figure 1: Map constituencies of groups
#####

## read in shape file
shp08 <- readShapePoly('01_data/03_maps/Constituencies.shp') # to open shape file

## export data table
df <- shp08@data
df$gis.id <- 1:nrow(df)

## make list of constits in experiments
sample <- subset(dat2, select = c(constit, province))
sample$harare <- ifelse(dat2$province=="Harare", 1, 0)
sample$constit <- gsub(' [0-9]$', "", sample$constit)
sample$constit <- toupper(sample$constit)


## change names of constits
sample$constit <- gsub(" - ", " ", sample$constit)
df$constit <- gsub(".", "", df$constit, fixed = T)
df$constit <- gsub("'", "", df$constit, fixed = T)
df$constit <- gsub("-", " ", df$constit, fixed = T)
sample$constit <- gsub('ZIBAGWE', 'CHIRUMANZU ZIBAGWE', sample$constit)
sample$constit <- gsub('MUREHWA NORTH', 'MUREWA NORTH', sample$constit)
sample$constit <- gsub('MPOPOMA', 'PELANDABA MPOPOMA', sample$constit)
sample$constit <- gsub('PLUMTREE', 'BEITBRIDGE EAST', sample$constit)
sample$constit <- gsub('WATERFALLS', 'HATFIELD', sample$constit)


## subset to unique
sample <- unique(subset(sample, select = c('constit', 'harare'), is.na(sample$harare)==F))


sample$constit[!(sample$constit %in% df$constit)]
df$constit[!(df$constit %in% sample$constit)]

## make sample indicator
df$sample <- ifelse(df$constit %in% sample$constit, 1, 0)

## make harare indicator
harare <- c("BUDIRIRO", "CHITUNGWIZA SOUTH", "CHITUNGWIZA NORTH", "EPWORTH", "DZIVARASEKWA",
            "GLEN NORAH", "GLENVIEW NORTH", "GLENVIEW SOUTH", "HARARE CENTRAL",
            "HARARE EAST", "HARARE NORTH", "HARARE SOUTH", "HARARE WEST", "HATFIELD",
            "HIGHFIELD EAST", "HIGHFIELD WEST", "KAMBUZUMA", "KUWADZANA",
            "KUWADZANA EAST", "MABVUKU TAFARA", "MBARE", "MOUNT PLEASANT",
            "MUFAKOSE", "SOUTHERTON", "ST MARYS", "SUNNINGDALE", "WARREN PARK",
            "ZENGEZA EAST", "ZENGEZA WEST")
df$harare <- ifelse(df$constit %in% harare, 1, 0)

## make bulawayo indicator
byo <- c("BULAWAYO CENTRAL", "BULAWAYO EAST", "BULAWAYO SOUTH", "EMAKHANDENI-ENTUMBANE",
         "LOBENGULA", "LUVEVE", "MAGWEGWE", "MAKOKOBA", "NKETA", "NKULUMANE",
         "PELANDABA-MPOPOMA", "PUMULA")
df$byo <- ifelse(df$constit %in% byo, 1, 0)

## put back in shape
df <- df[order(df$gis.id),]
shp08@data <- df

## add province borders to map
prov <- readShapePoly("01_data/03_maps/zim_provinces.shp")

## plot whole zimbabwe
plot_zim <- ggplot() +
  layer_spatial(shp08, aes(fill = factor(sample)), show.legend = FALSE) + 
  scale_fill_manual(values = alpha(c("white", "grey"))) +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank()) +
  layer_spatial(prov, lwd = 1.2, fill = NA)

## plot harare and bulawayo as insets

har08 <- subset(shp08, shp08$harare==1)

plot_har <- ggplot() +
  annotation_spatial(har08) +
  layer_spatial(har08, aes(fill = factor(sample)), show.legend = FALSE) + 
  scale_fill_manual(values = alpha(c("white", "grey")))+
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank())

byo08 <- subset(shp08, shp08$byo==1)

plot_bul <- ggplot() +
  annotation_spatial(byo08) +
  layer_spatial(byo08, aes(fill = factor(sample)), show.legend = FALSE) + 
  scale_fill_manual(values = alpha(c("white", "grey")))+
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank()) 


pdf('03_graphs/map_sample_all.pdf')
ggdraw() +
  draw_plot(plot_zim) +
  theme(plot.margin = unit(c(0,4.5,0,0), "cm")) + 
  draw_plot(plot_har, x = 1, y = 0.48, width = .3, height = .3) +
  draw_plot(plot_bul, x = 1, y = 0.15, width = .3, height = .3) +
  geom_rect(aes(xmin = 1, xmax = 1.3, ymin = 0.5, ymax = 0.8), fill = NA, col = 'black') +
  geom_rect(aes(xmin = 1, xmax = 1.3, ymin = 0.2, ymax = 0.45), fill = NA, col = 'black') +
  geom_text(aes(x = 1.15, y = 0.77, label = "Harare"), size = 4) +
  geom_text(aes(x = 1.15, y = 0.42, label = "Bulawayo"), size = 4)
dev.off()



#####
## Figure 3: Maps of violence by constituency
#####

sok <- read.csv('01_data/01_raw/sokwanele_violence.csv')

## change sok constit names
setnames(sok, 'consitname', 'constit')
sok$constit <- toupper(sok$constit)

sok$constit <- gsub("MT ", "MOUNT ", sok$constit)
sok$constit <- gsub("GLEN VIEW", "GLENVIEW", sok$constit)
sok$constit <- gsub(" - ", "-", sok$constit)
sok$constit <- gsub("MARAMBA-PFUNGWE", "MARAMBA PFUNGWE", sok$constit)
sok$constit <- gsub("CHIRUMANZU-ZIBAGWE", "CHIRUMANZU ZIBAGWE", sok$constit)
sok$constit <- gsub("ST MARYS", "ST. MARY'S", sok$constit)
sok$constit <- gsub("MABVUKU/TAFARA", "MABVUKU-TAFARA", sok$constit)
sok$constit <- gsub("GOKWE SENGWA", "GOKWE-SENGWA", sok$constit)
sok$constit <- gsub("GOKWE SESAME", "GOKWE-SESAME", sok$constit)


## subset to unique province-constit
constits <- unique(subset(sok, select = c('province', 'constit')))
dup <- duplicated(constits$constit)
constits <- subset(constits, dup==F)
sok$province=NULL
sok <- merge(sok, constits, by = 'constit')

## collapse sok to constit
viol_list <- c('assaultedwithhands', 'assaultedwithweapon',
               'abducted', 'falanga', 'submerged', 'strangled', 'burnt',
               'hung', 'tortured', 
               'theftdestrctionproperty') 
sok$events <- rowSums(sok[ ,viol_list])
sok$events <- ifelse(sok$events>0,1,0)
sok <- data.table(sok)
sok <- sok[,list('cio'=sum(cid, cio), 'zrp'=sum(zrpuniformed, zrpplainclothes, zrpriot),
                 'warvet' = sum(zanupfwarveteran), 'supporter' = sum(zanupfsupporter),
                 'youth' = sum(zanupfyouth), 'zna'=sum(zna),
                 'events' = sum(events)),
           by = c('constit', 'province')]


## read in shape file
zim <- readShapePoly('01_data/03_maps/Constituencies.shp')

## export data table
df <- zim@data
df$gis.id <- 1:nrow(df)

## merge w sok data
df$constit[!(df$constit %in% sok$constit)]
sok$constit[!(sok$constit %in% df$constit)]

df <- merge(df, sok, by = 'constit', all.x = T)
df <- df[order(df$gis.id),]

df$cio <- ifelse(is.na(df$cio)==T, 0, df$cio)
df$zrp <- ifelse(is.na(df$zrp)==T, 0, df$zrp)
df$warvet <- ifelse(is.na(df$warvet)==T, 0, df$warvet)
df$supporter <- ifelse(is.na(df$supporter)==T, 0, df$supporter)
df$youth <- ifelse(is.na(df$youth)==T, 0, df$youth)
df$zna <- ifelse(is.na(df$zna)==T, 0, df$zna)

## calculate proportion
df$all <- df$cio + df$zrp + df$warvet + df$supporter + df$youth + df$zna
df$cio.prop <- df$cio/df$all
df$zrp.prop <- df$zrp/df$all
df$warvet.prop <- df$warvet/df$all
df$supporter.prop <- df$supporter/df$all
df$youth.prop <- df$youth/df$all
df$zna.prop <- df$zna/df$all


## create harare and bulawayo dummies
df$hre <- ifelse(df$constit %in% harare, 1, 0)
df$byo <- ifelse(df$constit %in% byo, 1, 0)

## set violent events to 0 if missing
df$events[is.na(df$events)] <- 0

## put back in map
zim@data <- df


## set colors
brks <- quantile(zim$events, seq(0,1,0.1), na.rm=T)
zim$events_cut <- as.numeric(cut(zim$events, breaks = 10))
cols <- c('white', rev(gray.colors(15)[1:length(brks)]))

## create province layer
prov.layer <- list("sp.lines", prov, col = "black", lwd=1.2)


## plot whole zimbabwe
plot_zim <- ggplot() +
  layer_spatial(zim, aes(fill = events)) + 
  scale_fill_gradientn(colors = cols,
                       guide = guide_colourbar(title.position = 'top', hjust = 0.5, 
                                               barwidth = 15)) +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        legend.position = 'bottom',
        legend.text = element_text(size = 14),
        legend.title = element_text(size = 16)) +
  # scale_shape(guide = guide_legend(title.position = 'top')) +
  labs(fill = "Violent Events") +
  layer_spatial(prov, lwd = 1.1, fill = NA) 


## plot harare and bulawayo as insets

har08 <- subset(zim, zim$hre==1)

plot_har <- ggplot() +
  layer_spatial(har08, aes(fill = events_cut), show.legend = FALSE) + 
  scale_fill_gradientn(colors = cols[1:4]) +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank()) 

byo08 <- subset(zim, zim$byo==1)

plot_bul <- ggplot() +
  layer_spatial(byo08, aes(fill = events_cut), show.legend = FALSE) + 
  scale_fill_gradientn(colors = cols[1:2]) +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank()) 


pdf('03_graphs/viol_map.pdf')
ggdraw() +
  draw_plot(plot_zim) +
  theme(plot.margin = unit(c(0,4.5,0,0), "cm")) + 
  draw_plot(plot_har, x = 1, y = 0.48, width = .3, height = .3) +
  draw_plot(plot_bul, x = 1, y = 0.15, width = .3, height = .3) +
  geom_rect(aes(xmin = 1, xmax = 1.3, ymin = 0.5, ymax = 0.8), fill = NA, col = 'black') +
  geom_rect(aes(xmin = 1, xmax = 1.3, ymin = 0.2, ymax = 0.45), fill = NA, col = 'black') +
  geom_text(aes(x = 1.15, y = 0.77, label = "Harare"), size = 5) +
  geom_text(aes(x = 1.15, y = 0.42, label = "Bulawayo"), size = 5)
dev.off()



#####
## Figure 3: maps of poverty by constituency
#####

dhs <- read.csv('01_data/02_clean/dhs_clean.csv')

dhs$constit2 <- toupper(dhs$constit2)

dhs$constit2 <- gsub('EMAKHANDENI ENTUMBANE', 'EMAKHANDENI-ENTUMBANE', dhs$constit2)
dhs$constit2 <- gsub('GOKWE ', 'GOKWE-', dhs$constit2)
dhs$constit2 <- gsub('MHONDORO ', 'MHONDORO-', dhs$constit2)
dhs$constit2 <- gsub('ZVISHAVANE ', 'ZVISHAVANE-', dhs$constit2)
dhs$constit2 <- gsub('MUREHWA', 'MUREWA', dhs$constit2)
dhs$constit2 <- gsub('ST MARYS', "ST. MARY'S", dhs$constit2)
dhs$constit2 <- gsub('PELANDABA MPOPOMA', 'PELANDABA-MPOPOMA', dhs$constit2)
dhs$constit2 <- gsub('MABVUKU TAFARA', 'MABVUKU-TAFARA', dhs$constit2)

dhs$constit2[!(dhs$constit2 %in% df$constit)]
df$constit[!(df$constit %in% dhs$constit2)]

df <- merge(df, dhs, by.x = 'constit', by.y = 'constit2', all.x = T)
df$wgt_cuts <- as.numeric(cut(df$wgt_z/100, breaks = 10, include.lowest = T))
cols <- gray.colors(10)

df <- df[order(df$gis.id),]

zim@data <- df 


## plot whole zimbabwe
plot_zim <- ggplot() +
  layer_spatial(zim, aes(fill = wgt_z/100)) + 
  scale_fill_gradientn(colors = cols,
                       guide = guide_colourbar(title.position = 'top', hjust = 0.5, 
                                               barwidth = 15)) +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        legend.position = 'bottom',
        legend.text = element_text(size = 14),
        legend.title = element_text(size = 16)) +
  labs(fill = "Weight-for-height z-score") +
  layer_spatial(prov, lwd = 1.1, fill = NA)

## plot harare and bulawayo as insets

har08 <- subset(zim, zim$hre==1)

plot_har <- ggplot() +
  layer_spatial(har08, aes(fill = wgt_cuts), show.legend = FALSE) + 
  scale_fill_gradientn(colors = cols[1:9]) +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank()) 

byo08 <- subset(zim, zim$byo==1)

plot_bul <- ggplot() +
  layer_spatial(byo08, aes(fill = wgt_cuts), show.legend = FALSE) + 
  scale_fill_gradientn(colors = cols[c(4,5,7,10)]) +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank()) 


pdf('03_graphs/pov_map.pdf')
ggdraw() +
  draw_plot(plot_zim) +
  theme(plot.margin = unit(c(0,4.5,0,0), "cm")) + 
  draw_plot(plot_har, x = 1, y = 0.48, width = .3, height = .3) +
  draw_plot(plot_bul, x = 1, y = 0.15, width = .3, height = .3) +
  geom_rect(aes(xmin = 1, xmax = 1.3, ymin = 0.5, ymax = 0.8), fill = NA, col = 'black') +
  geom_rect(aes(xmin = 1, xmax = 1.3, ymin = 0.2, ymax = 0.45), fill = NA, col = 'black') +
  geom_text(aes(x = 1.15, y = 0.77, label = "Harare"), cex = 5) +
  geom_text(aes(x = 1.15, y = 0.42, label = "Bulawayo"), cex = 5)
dev.off()




